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ABSTRACT 

We use a 3-D relativistic SPH (Smoothed Particle Hydrodynamics) code to 
study mergers of black hole - neutron star (BH-NS) binary systems with low mass 
ratios, adopting q = M NS /M BH ~ 0.1 as a representative case. The outcome of 
such mergers depends sensitively on both the magnitude of the BH spin and its 
obliquity (i.e., the inclination of the binary orbit with respect to the equatorial 
plane of the BH). In particular, only systems with sufficiently high BH spin 
parameter a and sufficiently low orbital inclinations allow any NS matter to 
escape or to form a long-lived disk outside the BH horizon after disruption. 
Mergers of binaries with orbital inclinations above ~ 60° lead to complete prompt 
accretion of the entire NS by the BH, even for the case of an extreme Kerr BH. 
We find that the formation of a significant disk or torus of NS material around 
the BH always requires a near-maximal BH spin and a low initial inclination of 
the NS orbit just prior to merger. 

Subject headings: binaries: close - - black hole physics - - stars: neutron - 
relativity — gamma rays: bursts 
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Introduction and Motivation 



1.1. Double Compact Objects as Gravitational Wave Sources 



Over the past two decades, the modelling of double compact objects (DCOs) has at- 
tracted special interest among theorists, mainly because such systems are expected to be 
strong sources of gravitational waves (GWs). Their inspiral and merger GW signals cover a 
wide frequen cy band, from ~ 10~ 4 — 10 _1 Hz for supermassive BH binaries of ~ 10 4 — 10 7 M Q 
( lArunl 120061 ) all the way up to ~ 1000 Hz for mergers of NS-NS binaries, providing poten- 
tial sources both for ground-based interferometers (LIGO, VIRGO, etc.) and space-based 
detectors (LISA). The insp iral signals can provi de information on the spins and masses of 
the compact objects (e.g., iPoisson fc Willi Il995l ). Moreover the merger signals from BH 
NS and NS-NS binaries can carry information about t he NS internal struc t ure and the 



equation of state (EOS) o f matter at nuclear densities (IFaber fc Rasiol l2000l : iFaber et al. 



200ll ; IFaber & Rasiol 120021 ) . Note that for BH-NS mergers with fairly massive BHs, with 
Mbh ^ 100 Mq, where the NS is expected to plunge into the BH as a whole, not much 
information on the NS EOS will be carried by the GW signal. In contrast, the merger of 
a NS with a stellar-mass BH (Mbh ~ 10 Mq) makes it possible for the NS to be disrupted 
outside the BH's innermost stable circular orbit (ISCO), and this will potentially enrich the 
GW signal with a lot of information on the detailed behavior of the NS matter. Even though 
the GW signals from the inspiral of NS-NS binaries are accessible to ground-based interfer- 
ometers, covering the frequency range ~ 40 — 1000 Hz, the signals from their final merger s 
will probably be lost in the high-frequency noise level (jVallisneril |2000| ; IFaber et al.l 120021 ) . 
On the other hand, the GW merger signals for typical BH-NS binaries with stellar-mass BHs 
are expected to lie well within the sensitivity band of LIGO, at frequencies ~ 100 — 500 Hz. 



Although double NS binaries have been observed (IThorsett &: Chakrabartylll999l : lBurgay et al 



20031 ) . BH-NS and BH-BH binaries remain undetected. Moreover, the few observed NS-NS 
binaries (binary pulsars with NS companions and one double pulsar) are subject to consider- 
able selection effects. Therefore it is not currently possible to infer much empirically about 
the general properties of DCOs based on the observed sample. Theorists rely instead on 
binary evolution and population synthesis mo dels to make predictions about the formation , 
evolution and properties of such binaries (e.g.. lBelczynski et al.ll2007l ; iNelemans et al.ll200ll ). 
These models give estimates of merger rates for DCOs and the corresponding detection 
rates for the various GW interferometers. NS-NS binaries are expected to merge with rates 
~ 1 — 145 Myr -1 per MWEG (Milky Way Equi yalent Galaxy) and the equivalent rate fo r BH- 
NS binaries is ~ 0.07 - 5 Myr" 1 per MWEG kirn et al.lkood : iBelczvnski et al.ll2007h . For 



NS-NS binaries the detection rate estimates are ~ (0.4 — 60) x 10 3 yr 1 and ~ 2 — 330 yr 



for LIGO and Advanced LIGO, respectively (IKim et al.ll2006l ). while for BH-NS binaries the 
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2 x KTV -1 and 0.7 - 40 yr" 1 (IBelczvnski et al.l 120071 ). 



equivalent rates are ~3x 10~ 3 
Although the detection rates are quite low for the current LIGO stage, these predictions are 
very promising for Advanced LIGO. One should reme mber that, altho ugh the merger rates 
for NS-NS binaries have been empirically constrained (IKim et al.ll2006l ). no such constraints 
have been set for BH BH and BH-NS mergers. The lack of complete understanding of binary 
star evolution can lead to merger rate estimates for these systems which vary significantly, 
depending on the exact physical assumptions adopted in the various binary evolution codes. 
Current effort is focusing on decreasing these uncertai nties and setting more solid constr aints 
on the merger rates of BH-NS and BH-BH binaries (jO'Shaughnessy et al.ll2005l 120061 ). 



Along with the binary evolution studies that provide merger rates of compact binaries, 
general relativistic calculations of binary mergers try to guide the search for GW signals 
by predicting the exact shape of the signals and generating GW search te mplates (e.g., 
Buonanno etaDl2007i : baker et al.lboOTl : I Abbott et aklhood : lApostolatoslll995h . For BH-NS 
and BH-BH binaries it is expected that both the BH spins and their possible misalignment 
angle with respect to the orbital an gular momentum will affect significantly the s hape of the 
GW signals and their detectability (1 Apostolatosl 1 1 9951 : iGrandclement et al.ll2003l ). 



1.2. Connection to Gamma Ray Bursts 

Another interesting aspect of BH-NS and NS-NS binaries is their possible connection 
to the observed short gamma-ray bursts (GRBs). This hypothesis has gained widespread 
support over the last few years, both because of the rapid progress in theoretical modelling 



and from the recent Swift observations of short GRBs (see lNakarll2007l . for a recent review). 



GRBs are classified into two duration classes, separated at ~ 2 s (Kouveliotou et al. 
1993). Long bursts are found to be predominantly in active star-forming regions. It is now 
believed that long bursts are produced when a massive star reaches the end of its life, its 
core collapsing to form a BH and, in the process, ejecting an ultra-relativistic outflow (e.g., 
Woosley & Bloom 2006). The standard collapsar model predicts that a broad-lined and 
luminous Type I-c core collapse supernova (SN) accompanies long bursts (MacFadyen & 
Woosley 1999). This association has been confirmed in observations of several nearby GRBs 
(e.g., Galama et al. 1998; Hjorth et al. 2003; Pian et al. 2006; Gehrels et al. 2006). 

Until recently, afterglow of short bursts have been extremely elusive. This situation 
changed dramatically in 2005. Swift and HETE-2 detected X-ray afterglows from short 
bursts (Gehrels et al. 2005; Villasenor et al. 2005). This has led to the identification of 
host galaxies and to redshift measurements. More than 10 short burst afterglows have been 
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detected so far, and distinctive features have emerged. While long bursts occur only in star 
forming spiral galaxies, short bursts appear also in elliptical galaxies, which are dominated 
by an old stellar population. The low level of star formation makes it unlikely that the 
bursts originated in a SN explosion. Even though a short burst, GRB 050709, was seen in 
a galaxy with current star formation, optical observations ruled out a SN association (Fox 
et al. 2005). The isotropic energy for short bursts is 2-3 orders of magnitude lower than 
that for long bursts E iso ~ 10 52 ~ 54 erg (Barthelmy et al. 2005). These results suggest that 
compact stellar mergers are the progenitors of short bursts. 

The similarity of X-ray afterglow light curves of long and short bursts indicates that 
afterglows of both classes can be described by the same paradigm, despite the differences 
in the progenitors. This view is supported by the fact that the decay rate of short burst 
afterglows is the value expected from the standard fireball model (e.g., Piran 2004; Meszaros 
2006), and that at least in two short bursts (GRB 050709 and GRB 01221 A) there is evidence 
for a jet break (Fox et al. 2005; Soderberg et al. 2006; Burrows et al. 2006). In the standard 
afterglow model, these breaks are interpreted as a signature of collimation of a fireball into 
a jet with an opening angle 9 ~ 6 — 12 degrees and imply a beaming-corrected energy of 
E ~ (0.5 — 3) x 10 49 erg, much less than that of long bursts, which have E ~ 10 51 erg (Frail 
et al. 2001). The lower energy implies that the mass of the debris torus formed during the 
merger could be smaller than that of the torus formed in the collapse of the core of massive 
stars. 

Combined with the lack of a jet break in GRB 050724, which gives lower limits of 9 > 25 
degrees and E > 4 x 10 49 erg (Grupe et al. 2006), the current small sample indicates that 
the outflow of short bursts is less strongly collimated than most previously reported long 
GRBs with the median value 9 ~ 5 degrees (Frail et al. 2001; however see Monfardini et al. 
2006). The wider jet angle is consistent with a merger progenitor scenario (e.g., Meszaros, 
Rees & Wijers 1999), since there is no extended massive stellar envelope (as in long GRBs) 
that serves to naturally collimate the outflow. Many more bright short bursts will be needed 
to improve the jet break statistics substantially. 

One of the unexpected results from Swift is that early X-ray afterglows of long bursts 
show a canonical behavior, where light curves include three components: (1) a steep decay 
component, (2) a shallow decay component and (3) a "normal" decay component. On top 
of this canonical behavior, many events have superimposed X-ray flares (e.g., Zhang et al. 
2006). The X-ray afterglow of a short burst, GRB 050724, associated with an elliptical host 
galaxy, also resembles the canonical light curve, and it suggests a long-lasting engine. A 
flare at ~ 100 s in the X-ray light curve decays too sharply to be interpreted as the afterglow 
emission from a forward shock, but is consistent with the high latitude emission from a 
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fireball (Barthelmy et al. 2005; Kumar & Panaitescu 2000). This is appropriate for the late 
internal shock scenario as invoked to interpret X-ray flares in long GRBs. This interpretation 
requires that the central engine remains active up to at least ~ 100 s, and challenges simple 
merger models, because the predicted typical time scale for energy release is much shorter. 

Another interesting scenario based on BH-NS and NS-NS mergers and connected to 
GRBs is the production of r-process elements (heavy nuclei with A > 90 — 100) through the 
nuclear physics of decompressed NS matter in the ejecta produced by the merger. It is still 
not clear today what is the astrophysical site that can provide the appropriate conditions for 
r-process nucleosynthe sis to take place, although the conditions themselves can be estimated 
( IJaikumar et al.ll2006l ). The possible ejection of extremely neutron-rich (Y e ~ 0.1) material 
from NS disruption s in compact binary mergers is b elieved to be a promising source for 



r-process elements ( Lattimer Sz Schramm 1974. 1976h . As s uch mergers are expected to 



happen in the outskirts of galaxies (IPerna fc Belczynskill2002l ). it is possible that the high- 
velocity electa will also enrich the intergalactic medium with high mass r-process elements 



(IRosswog 



2005|). 



1.3. Results from Previous Studies 



In recent years, various groups have performed 3-D hydrodynamic simulations of BH-NS 
mergers, providing very interesting — if somewhat preliminary — results. In parallel with this 
work there has been an ongoing effort to improve both the computational techniques and 
the overall accuracy of the merger simulations. 

Newtonian studies of BH-NS mergers, w here the BH is simply represented by a point 



mass, were first used by different groups (e.g., iLee fc Kluzniaklll999bl ; IRosswog et ajj 120041 ) 
to cover a fairly wide range of mass ratios while attempting to determine the merger's depen- 
dence on the NS EOS. The general picture that comes out of these Newtonian calculations 
is a clear connection between mass ratio and NS EOS, and the ability to form a massive 
enough accretion disk around the BH at the end of the merger. As an overall trend, for 
higher mass ratios (with q ~ 0.7 — 0.28) the survival of NS material and formation of a disk 
around a Schwarzschild BH seems possible for soft EOS. For these high mass ratio merg- 
ers the tidal disruption radius lies outside the ISCO, which allows for more NS material to 
settle into a disk after the tidal disruption. However, those studies also suggested that, for 
a stiffer EOS, the core of the NS always survives the merger, inhibiting the formation of a 
massive disk around th e BH and leading instead to a period of quas i-stable mass transfer 
from the NS to the BH Lee fe KluzniakHl999al : iKluzniak feTeelll998h . A stiffer EOS could, 



in that case, prevent the total disruption of the NS, leading to the eventual formation of 



- 6- 



a "mi ni-NS" along with th e disk or even lead to multiple disruptions of the surviving NS 



core. 



Rosswog et al.l (120041 ) suggested that the survival of a mini-NS for the case of a stiff 
polytropic EOS (r = 3) is connected to the difficulty of forming a massive disk in their 
calculations, arguing that the survival of an orbiting NS core acts as a storage mechanism 
that prevents further inflow of material towards the BH. 

In their full-GR tr e atmen t of BH-NS mergers for a non-spinning BH of arbitrary mass, 



Shibata fc Uryul (120071 . 120061 ) concluded that the disruption of a NS by a low-mass BH 
(M = 3.2 — 4M ) can lead to the formation of a low-mass disk (of mass ~ 0.1 M & ) around 
the BH, which could potentially power a short GRB. No formation of more massive disks 
was ever observed in their simulations, indicating that systems with ~ M & disks around a 
BH cannot be formed through BH-NS mergers with non-spinning BHs. Survival of a NS 
core, or quasi-stable mass transfer, as predicted by Newtonian calculations, are never seen 
in these full-GR calculations. 



Also limited to Schwarzschild BHs, the fully relativistic calculations of iFaber et al. 



( ]2006bl ) focused on BH-NS mergers with a much more extreme mass ratio q = 0.1. For 
such a mass ratio and a non-spinning BH, the tidal disruption limit of a NS with canonical 
mass and radiu s (i.e., compactness ratio C = M^s/Rns — 0.2) is reached inside the ISCO. 



For this reason, IFaber et al.l (l2006bl ) considered artificially undercompact models for the NS 
(with compactness as small as C ~ 0.04), in order to study cases where the disruption of 
the NS takes place outside the ISCO. They suggested that their study for these undecom- 
pact NSs can serve as an analogue for binaries with lower-mass BHs and more compact NSs 
(which their numerical code could not handle directly), where the tidal radius is located well 
outside the ISCO. They found that only a small fraction (~ 5 — 7%) of the NS mass becomes 
unbound and escapes from the system and that, while most of the infalling mass is accreted 
promptly by the BH, part of it (~ 25%) remains bound outside the horizon, forming a disk. 

Motivated by all these recent observational and theoretical developments, we have em- 
barked on a new numerical investigation of the merger process for BH-NS binaries using a 
3-D relativistic SPH code developed specifically for the study of stellar disruptions by BHs. 
In this paper, the first of a series, we present our methods and numerical code, as well as the 
results from a first set of preliminary calculations aimed at exploring broadly the parameter 
space of these mergers. Our paper is organized as follows: In §2 we describe the SPH code 
used for our simulations; we develop some analytic considerations regarding the metric used 
by the code; and we also discuss the various test calculations that we have performed. In 
§3 we present results from our simulations of equatorial mergers (§3.1) and inclined mergers 
(§3.2), including a discussion of how to set up initial conditions for the inclined case. The 
GW signals extracted from these simulations are presented and discussed in §4. Finally, a 
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summary and conclusions are given in §5. 

2. Methods and Tests 

2.1. Critical Radii 

Fig. [H shows why the final merger of a BH-NS binary (for a typical stellar-mass BH) 
is interesting but also particularly difficult to compute: the tidal (Roche) limiiQ is typically 
right around the 1SCO and the BH horizon. On the one hand, this implies that careful, fully 
relativistic calculations are needed. On the other hand, it also means that the fluid behavior 
and the GW signals could depend sensitively (and carry rich information) on both the masses 
and spins of the compact objects, and on the NS EOS. How much information is carried about 
the fluid depends on where exactly the tidal disruption of the NS occurs: for a sufficiently 
massive BH, the horizon will always be encountered well outside the tidal limit, in which 
case the NS behavior remains point-like throughout the merger and disruption will never be 
observed. In the opposite case where the tidal limit resides outside the horizon, the GW 
signal corresponding to the NS disruption could be detected by ground-based interferometers. 

For those cases where disruption occurs outside the BH's horizon, the final outcome of 
the merger depends strongly on the relative positions of the tidal radius R t (i.e., the point 
where disruption takes place) to the ISCO of the BH. Fig. [1] gives an overview of these 
critical radii and a first idea of what one could expect the outcome to be for mergers with 
various mass ratios. An interesting fact that we notice first is that for a given BH mass the 
relative position of Rt to the ISCO changes with the BH's angular momentum: as the BH 
spin increases, it drags the ISCO closer to the BH. For high enough mass ratios (low BH 
mass of the order of a few M ) the tidal radius is always encountered outside the ISCO, 
even for non-rotating BHs. For a 10 M & BH, the ISCO and R t coincide for a Schwarzschild 
BH, but the ISCO moves inside R t for a Kerr (spinning) BH. Finally, for even higher mass 
BHs of ~ 15 M Q the situation becomes more interesting as the ISCO lies well outside the 
tidal radius for a non-spinning BH and for Kerr BHs with low spin, but it starts migrating 
inside R t as the BH's angular momentum increases. For fast spinning BHs the tidal radius 
is encountered well outside the ISCO. Of course for very massive BHs, R t is not only inside 
the ISCO but inside the horizon as well, even for the case of an extremal Kerr BH (this 
happens for M > 100 M ). 



1 see, e.g., 
radii. 



Lai et al. 



( 19941 ) and IWiggins fc Lail ( 200(1 ) for analytical calculations and discussion of tidal 
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These simple results are important because disruption outside the ISCO could lead to 
the formation of a disk of NS debris outside the BH's horizon. If the disruption is to happen 
inside the ISCO, no such feature is expected. Knowing the relative position of ISCO and 
R t for various binary mass ratios therefore gives us a first insight on what to expect as 
the qualitative outcome of a merger, and in which ranges of mass ratios we should look for 
certain outcomes. 



2.2. Analytic Considerations 



The SPH code makes use of the Kerr-Schild (K-S) form of the Kerr metric. In this 
section we summarize the reasons for choosing to use the K-S metric (i.e., its advantages 
over the Kerr metric in Boyer-Lindquist [B-L] coordinates) and how quantities such as the 
BH horizons and the angular velocity of equatorial circular orbits around the BH translate 
from one coordinate system to the othe r. The reader can fin d extensive disc u ssions of th e 
B-L and K-S coordinate systems in, e.g.. IChandrasekharl (119831 ) ; |Poissonl (12004 ) ; iKerrl (119631 ). 
Here we present a brief overview with emphasis on some useful aspects of the K-S metric 
that relate directly to our calculations and to the presentation of our numerical results. 



The Kerr solution in B-L coordinates is given by the familiar expressions ( IBoyer-Lindquist 



19671 ) 



ds 2 



9 Mr y 

;i _ ^)dt 2 + fydr 2 + Ed9 2 



. 2 2 2Mra 2 . , „, 
+ (r + a H sm & sm' 



, „ , , , AMra sin 2 9 , , , 
- a1 - 2 do dl 



where 



A = r 2 - 2Mr + a 2 
E = r 2 + a 2 cos 2 9. 



(2) 
(3) 



This form has only one off-diagonal term and is therefore far more convenient to use 
than the K-S form of the metrical Yet, it carries some extra coordinate singularities which 
correspond to the roots of A : r± = M ± y/ M 2 — a 2 (r + and r_ are the future and past 
horizons respectively, with r + only being an event horizon). It is useful to observe that: (a) 



2 Note that this form of the Kerr metric reduces to the well known Schwarzschild solution in the limit 
a = 0: ds 2 = -(1 - ^)dt 2 + (1 - —)~ 1 dr 2 + r 2 d9 2 + r 2 sin 2 ^ 2 



for a/M = 1 the two coordinate singularities (horizons) coincide, (b) for a/M = there 
is only one horizon at r^/M = 2 and the curvature singularity at r/M = 0. The obvious 
advantage of casting the Kerr metric into its K-S form is that one avoids the coordinate 
singularities at the horizon present in B-L coordinates. 

The coordi nate singularit y at the horizon present in the B-L form of the metric has the 



following effect (jPoissonl 12004 ): although it takes a finite proper time for a particle to cross 
the event horizon, it takes infinite coordinate time t to do so. Moreover, since the angular 
velocity dcj)/dt tends to a finite limit at the horizon, has also to increase an infinite amount 
: <p — > oo (<f), like t, is not a "good" coordinate at the horizon). What this practically means 
for our code is that it prevents us from extending our calculations all the way to the horizon 
of the BH. In that case, one needs instead to place an absorbing boundary outside the BH's 
horizon. Fortunately there is a way to overcome this problem. In order to extend the Kerr 
metric beyond the horizon, another coordinate system needs to be adopted. Keeping in 
mind that the horizons are null surfaces, it makes intuitively sense to construct the new 
coordinates in terms of null geodesies. 

The null geodesies in the Kerr spacetime are given by the tangent vectors 

dt r 2 + a 2 



dr A 



-E (4) 



dr 



where E is the specific energy. 
The real null vector I reads: 



l<* = -^(r 2 + a 2 ,±A,0,a). 
Setting E = 1 (and using A for the affine parameter) 



d\ A ' dX ' [ ) 
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By choosing the positive sign for dr/dX we obtain an outgoing congruence with the 
tangent vector field defined by 



r 2 + a 2 . a 



l a $o l = —^- tft + tfr + ^V (11) 



The new variables u and can be introduced in the place of t and 4>: 

2 2 

du = dt- T ^ Q dr (12) 
d^) = d<p- ^dr. (13) 

The null geodesic now becomes 

l a = (0,1, 0,0) (14) 

and the metric takes the form 

, o 2Mf, , 2 o , 4aMf sin 2 6 1 , ~ , ,.,„ - 
as = — (1 ^^)du + Latf — zdudr — au a(p + 2a sin Udrdcp 

, 2 o 2Mfa 2 sin 2 6' . 2 „,72 
+ (r 2 + a 2 + ) sin 2 6d<p 2 (15) 

Zj 

where f is defined bj§| 

- (p 2 - a 2 )r~ 2 - aV = (16) 

witlfl 

p 2 = x 2 + 2/ 2 + z 2 . (17) 

By instead choosing the negative sign for dr/dX we can obtain an ingoing congruence 
with the tangent vector field given by 

l' a #a = ^^-#t-#r + fa. (18) 



5 a more geometrically insightful representation of Eq. (|16[) would be 



x 2 + y 2 i z 2 _^ 



'Note that as a — > , f — » p. 
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The new variables to be introduced here are 



r 2 + a 2 



dv = dt-\ dr (19) 

dfi = d(p+ ^dr. (20) 
In this coordinate system the null geodesic simplifies to 

r = (0,-1,0,0) (21) 

and the metric becomes 

2 Mr An Mr sin 2 R 

ds 2 = -(\- —)dv 2 + Y>dR 2 + 2dvdr - - dvdfi - la sin 2 Rdfd$ 

, 9 9 2Mra 2 sin 2 R . 9 „ , ~ l0 , . 

+ (f 2 + a 2 + ) sin 2 ^d0' 2 . (22) 

What is the use of those two different coordinate sets? The coordinates (v,r,R,(j)') are 
well behaved on the future horizon, yet they are singular on the past horizon, where (u,r,R,<j)) 
are now well behaved. Therefore the first set of coordinates is used in order to regularize 
the past horizon, whereas the second one is used to regularize the future horizon. If one 
wants to avoid both horizons of a Kerr BH, both patches need to be used in order to cover 
the entire spacetime around the BH. Then there is just the curvature singularity of the Kerr 
metric left, which occurs at 

S = f 2 + a 2 cos 2 9 = 0, (23) 



implying 



x 2 + y 2 = a 2 . (24) 



The curvature singularity is not a point but rather a ring of radius a and exists only in the 
equatorial plane. It is always found inside the past horizon, although as a/M decreases the 
curvature singularity and the past horizon approach each other (Fig. [2]) . 



2.3. SPH Code 

We employ a 3-D implementation of the Lagrangian SPH technique. Our GRSPH 
(General Relativistic Smooth Particle Hydrodynamics) code is based on the work by Laguna 
et al. (1993a) in which the general relativistic hydrodynamic equations were rewritten in 
a Lagrangian form similar to their counterparts for non-relativistic fluids with Newtonian 
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self-gravity. The GRSPH code is restricted to fixed curved spacetimes; in particular, the 
spacetime of a rotating BH is assumed here. The use of a fixed background in our simulations 
is justified for low mass ratios (q ~ 0.1 in this paper), i.e., when the mass of the BH is 
substantially larger than the NS mass. 

Originally the BH metric used in the code was in terms of Boyer-Lindquist coordinates. 
For the present work, we use Kerr-Schild coordinates. The main advantage of these coordi- 
nates is their regularity across the BH horizon, thus allowing SPH particles to freely cross 
it. An important aspect of SPH in curved spacetimes is the handling of the local volume- 
averaging required by the smoothing of the equations. The smoothing volumes involved are 
in general not small enough to ignore curvature effects, especially in the neighborhood of the 
BH. Our implementation correctly accounts for these effects. The most important compu- 
tational aspect of the GRSPH code is the particle neighbor findin g algorithm. The GRSPH 
code uses an oct-tree data structure as basis for finding neighbors (IWarren fc Salmonlll99ll ). 
The oct-tree neighbor- searching p art of GRSPH has als o been successfully used for iV-body 
large-scale structure simulations (IHeitmann et al.ll2005l ) and is used as the foundation of a 
different SPH cod e with radiation transport used for supernova core collapse calculations 
rtFrver et alJhood ). The GRSPH code scales as 0{N log N) with N the total number of SPH 
particles. The code was calibrated with three one-dimensional benchmarks (Laguna et al. 
1993): (1) relativistic shock tubes, (2) dust infall onto a BH, and (3) Bondi collapse (see 
Laguna et al. 1993a for details). The code has been successfully applied in several previ- 
ous studies of the tidal disruption of ordinary (main-sequence) stars by a supermassive BH 
(Laguna et al. 1993b; Kobayashi et al. 2004; Bogdanovic et al. 2004). 



m 



Artificial viscosity i s imp lemented in our cod e by following the description presented 



Laguna et al. 


( 


1993a 


) and 


Laeuna et al. 


f 


1994) 



account for the possible prese nce of shocks and is introduced as a viscous pressure term 
added to the SPH equations (ILaguna et al.l Il993al ). As we do not expect any shocks to 
appear in our simulations except perhaps at the point where the stream of NS material 
accreting onto the BH self-intersects, we keep artificial viscosity suppressed by lowering the 
parameters of its two terms (the bulk viscosity term and the von Neumann-Ric htmyer (1950) 
viscosi ty term) to the value 0.2 instead of the more common values ~ 1 (see 
1993al . for details and tests). 



Laguna et al. 



Radiation reaction is implemen ted in the GRSP 
proximate prescription described in iLee fc Kluzniakl ( 



i code by following the simple, ap- 
1999bJ), §2. We use the quadrupole 



formula for th e rate of energy radiated (Eq. 4 in lLee fc Kluzniaklll999bl ) to derive a damping 
force (Eq. 6 in lLee fc Kluzniakl Il999bl ) for each SPH particle. The formula has been slightly 
modified as in our code we deal with du/dt where u is the 4-velocity. We have also added a 
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parameter to increase the radiation reaction force in order to accelerate the initial inspiral 
since we do not want to spend too much computational time on this initial phase. We have 
checked that, once the merger starts, the precise value of this parameter does not affect our 
results significantly. Throughout our calculations and for all the results presented here we 
use geometrized units with G = c = 1. 



2.4. Test Calculations 

2.4-1- Testing the K-S Metric 

The Kerr-Schild coordinate system used in the SPH code is only avoiding the outer 
(future) horizon, since that serves the purpose of the SPH particles getting as close to the 
BH's horizon as possible, without the code crashing or becoming problematic. Yet, in the 
extremal Kerr case (a/M = 1), the two horizons coincide (Fig. |2j) at f = \/2M, and as a 
result, we see the orbits of the SPH particles being trapped at f = \/2M. To check that 
it is the inner horizon where our transformation is singular, we set up the following simple 
test: we construct a test particle geodesic integrator that makes use of the K-S metric and 
we start with an equatorial orbit for a/M = 1 which, starting from a finite distance from the 
BH (having v z = v y = y = z = 0), ends being trapped at f = \/2M. We then keep following 
the same orbit as we reduce the value of a (Fig. [3]) . The result is that the particle ends being 
trapped at the inner horizon (whichever that is for the specific value of a). The trapping 
is illusionary: it reveals the singularity of our coordinate transformation. We thus avoid 
using a/M = 1 in the simulations presented in this paper and instead we set a/M = 0.99 
for the case of an extremal Kerr B.H. By doing so, the future horizon moves outwards at 
f_|_ = 1.51067M and the past horizon moves inwards at f_ = 1.31067M and are therefore 
distinguishable. 

Another quantity that is going to be affected by the coordinate transformation (i.e., 
from B-L to K-S) is the angular velocity for a circular equatorial orbit around the BH. 
In B-L coordinates is given by 



n + ~ ± r 3/2 ±aM l/2 ( 25 ) 

where r 2 = x 2 + y 2 + z 2 and the upper (lower) sign corresponds to co-rotating (counter- 
rotating) orbits. In K-S coordinates Eq. ( 125]) holds when r is replaced by f, as defined by 
Eq. ( |T6i) (see Appendix C for an analytic calculation of in K-S coordinates). In order 
to check this result numerically, we set up the following test: we use again the geodesic 
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integrator and we place a test particle at some distance ro from the BH and give it the 
angular velocity that corresponds to an equatorial circular orbit at the specific r as required 
by Eqs. (|25p and f[T6"j) . As expected, the test particle remains in this fixed circular orbit 
(with dr/dr = 0) as we integrate for a large number of orbital periods^. 



2.4-2. Stable Binary Orbits 

In a first test, we considered a white dwarf (WD) with M* = 0.6 M & and R* = 1.3 x 
1CT 2 -R Q ~ 9 x 10 8 cm orbiting around a Schwarzschild BH with M = 2 x 1O 5 M . Here 
5000 particles are used to represent the WD with a T = 5/3 polytropic EOS. Since the mass 
ratio is extreme, our approximation (moving the fluid on a fixed background metric) should 
be extremely accurate. With these parameters, the tidal radius Rt ~ 6.3 x 10 10 cm and the 
horizon scale = 2M ~ 5.9 x 10 10 cm are comparable, and much larger than the WD radius 
R*. If R t ^> Th, the point particle approximation should break down for orbits near the 
ISCO because the WD will get disrupted at radii well outside the ISCO. On the other hand, 
when R t <C rh, the sound crossing time of the WD is much shorter than the orbital period, 
and it is numerically expensive. Thus, we have chosen the parameters satisfying R t ~ to 
test the code for relativistic orbits. With these parameters we found that we can maintain 
a circular orbit at r = 8M to within \Ar\/r < 10~ 3 over one full orbital period (Fig. H|). 

In a second, similar test we considered a NS with M^s = 1-4 M and Rns = 13.4 km 
= 1.93 x lO" 5 ^ (represented by 10000 SPH particles and with a Y = 2 polytropic EOS) 
orbiting around a Schwarzschild BH with M = 10 M . We found that we can maintain a 
NS orbiting at r = 20 M stably and without any noticeable numerical dissipation for more 
than 20 orbital periods. 



2.5. Initial Conditions for BH— NS Binaries 



We set up initial conditions for BH-NS binaries near the Roche limit using the SPH 
code and a r elaxation technique sim ilar to those used for previous SPH studies of close 
binaries (e.g.. lRasio fc Shapiro! 1 1 99 5l ). First we construct hydrostatic equilibrium NS models 
for a simple gamma-law EOS by solving the Lane-Emden equation. When the NS with this 
hydrostatic profile is placed in orbit near a BH, spurious motions could result as the fluid 



5 If we instead use Eq. (|25| with r 2 = x 2 + y 2 + z 2 , the particle follows an oscillating orbit (dr/dr ^ 0) 
around its initial position, indicating that the formula used for needs modification. 
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responds dynamically to the sudden appearance of a strong tidal force. Instead, the initial 
conditions for our dynamical calculations are obtained by relaxing the NS in the presence of 
a BH in the co- rotating frame of the binary. For synchronized configurations (assumed here), 
the relaxation is done by adding an artificial friction term to the Euler equation of motion 
in the co-rotating frame. This forces the system to relax to a minimum-energy state. We 
numerically determine the angular velocity Q corresponding to a circular orbit at a given r as 
part of the relaxation process. The advantage of using SPH itself for setting up equilibrium 
solution is that the dynamical stability of these solutions can then be tested immediately by 
using them as initial conditions for dynamical calculations. 



3. First Results 

3.1. Equatorial BH— NS Mergers for Spinning and Non-spinning BH 

For all the results presented in this section we ran simulations using N = 10 4 SPH 
particles to represent a NS with a T = 2 polytropic EOS, except for one case, for which we 
have performed a crude convergence test by repeating the run with N = 10 5 particles (Run 
El of Tabled]). The number of neighboring particles used in the SPH code is = 80 and 
Nn = 140 for the low and high resolution runs, respectively. The NS mass is M^s — 1.4M 
and its radius is Rns — 15 km. The BH has a mass M = 15M . All SPH particles are of 
equal mass {m p = M NS /N). The NS is initially placed at a distance r = 8M and for cases 
with a spinning BH (a/M ^ 0) the NS is co-rotating in the equatorial plane. The angular 
momentum of the BH is left as a free parameter. 

In this first study we would like to get an idea of the possible outcomes of a BH-NS 
merger: what percentage — if any — of the NS mass survives the merger; what are the mor- 
phological features of the merger (e.g., creation of an accretion disk, spread of the accreting 
material), and how the BH angular momentum affects these features. For the two extreme 
cases a/M = and a/M = 0.99 (we avoid using a/M = 1 for the reasons mentioned in §2.2), 
we observe completely different behaviors. 

Fig. shows the outcome of the merger for a non-spinning BH: the NS, after being 
completely disrupted and following an inspiraling orbit, disappears entirely into the BH's 
horizon in a time t/M ~ 180 after the beginning of the simulation (where t is the coordinate 
time for an observer at inifinity). 

The result of our simulation for the case of the extremal Kerr BH (a/M = 0.99) is shown 
in Figs. [6] and [7J The NS again gets completely disrupted and falls toward the co-rotating 
BH following an inspiraling orbit while an outwards expanding tail of the disrupted NS's 
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material is forming at the same time. The NS fluid starts disappearing into the BH's horizon 
at t/M ~ 400 with an initially high infall rate which finally diminishes to an almost zero 
level. At this point {t/M = 490) about 33% of the initial NS mass resides outside the BH's 
horizon (Fig.|SJ). We notice that the disruption takes place well outside the ISCO, in contrast 
to the previous case of a Schwarzschild BH. 

We run the same exact type of simulations for different values of the BH's angular 
momentum and compare the results. Fig. [H] shows the fraction of the initial NS mass that 
survives the merger for five different values of a/M (0.75 < a/M < 0.99). As the angular mo- 
mentum of the BH decreases, so does the fraction of material that survives, residing outside 
the BH's horizon. For the case of a/M = 0.1 (not included in the latter figure for practical 
reasons), the situation is almost the exact same as for the Schwarzschild (non-spinning) BH: 
the NS disappears completely into the horizon (Table [T] summarizes our results). 

For the simulations with a/M spanning the range 0.1 < a/M < 0.75, namely for 
a/M = 0.2, a/M = 0.5 and a/M = 0.6, we observe that the infall starts earlier as a/M 
decreases (consistent with the behavior of the high a/M mergers). No surviving material 
exists for those mergers with low spin for the BH (a/M < 0.7). 

The morphological features of the mergers for the non-maximally spinning BH are simi- 
lar to the extremal case, with the difference that, as a/M decreases the outwards expanding 
tail tends to spread less. For the very low a/M cases no such tail is forming. 

By the end of the SPH simulations, the mass fraction that resides outside the BH's future 
horizon seems to have reached a stable value. Since for a/M > 0.9 there is a substantial 
fraction of the NS mass surviving the merger, one would like to investigate what the final 
fate of this material might be, i.e., what percentage of it — if an — will escape or stay bound, 
forming a stable disk around the BH. To answer those questions we first calculate the rest 
energy (energy as measured by an observer at infinity) of each SPH particle throughout 
the whole simulation for runs E1-E5. By knowing a particle's energy we can determine 
whether it is bound or unbound (e 2 > 1 for unbound and e 2 < 1 for bouncj^l, where e = 
— (<7ttW* + gt^u^ + gt r u r ) is the energy per unit rest mass for an observer at infinity). Fig. [9] 
shows the variations of bound and unbound mass for each of our simulations. 

For the two runs with the highest a/M values (runs El and E2 from Table 1) we are 
able to resolve a nonzero fraction of bound material (2 — 3% of the NS mass). For runs 
E3, E4 and E5 the percentage of bound material drops to unresolvable values (recall that 



3 for an analytical discussion and a proof of that see IWilkinsI (|l97lh and Ichandrasekhail (|l983l chap. 3, 



jl9) 
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our mass resolution is m p = M^vs/lO 4 ). As shown in Fig. [9] for the runs El and E2, both 
the unbound and bound NS material stabilize at a certain value and persist there for a 
long time. That leads us to believe that the bound fraction of the surviving mass stays 
around the BH forming a stable disk. In order to check whether our definition of bound 
and unbound material is correct, we perform the following test-runs: we evolve the results of 
the SPH simulations using again the geodesic integrator. For every run (different a/M) we 
set as initial conditions for the geodesic integrator the last output file from the equivalent 
SPH simulation. We make sure to run the geodesic integrator for a sufficiently long time, 
as it is much faster than the SPH code. The results of those tests confirm our expectations. 
Namely, the material that we recognized as unbound escapes completely following parabolic 
trajectories, whereas the particles that were energetically determined to be bound (for runs 
El and E2 only) remain bound around the BH (and outside the BH's horizon) following 
equatorial precessing orbits, with maximum apastron ~ 30M, minimum periastron of about 
~ 2 — 3M and a very small dispersion in the direction of the BH spin axis ~ 0.1 — 0.2M (see 
Fig. HO]). 

Two important results are seen in Table [U no significant fraction of bound material 
survives for mergers with a/M < 0.95 and therefore no formation of a stable disk is observed; 
for the cases with a/M < 0.7 the merger is completely catastrophic for the NS, with no 
material surviving. 

Finally, in the upper two plots of Fig. M we present a comparison between a lower 
resolution run (10 4 particles) and a run with higher resolution (10 5 particles) for the case 
of a maximally spinning BH. The left panel corresponds to the low-resolution run and the 
right panel to the high-resolution run. The agreement between the two runs is very good and 
this justifies our use of the lower resolution for all runs presented in this paper. The higher- 
resolution run took 248 CPU hours on two dual-core AMD Opteron 2.8 GHz processors, 
while the lower-resolution run took only 26 CPU hours. For both cases, the simulation 
resulted in the ejection of about 30% of the NS mass (unbound material) along with the 
survival of ~ 2 — 3% of the initial NS mass as bound material. 



3.2. Non-equatorial BH— NS Mergers 



Many studies have suggested t hat a significant NS birth kick is requ ired for the formation 



of coalescing BH-NS binaries (e.g.. lKalogerall2000l ; iLipunov et al.lll997l ). Any misalignments 



between the axis of the BH and the NS progenitor are expected to be canceled during the 
evolution of the binary prior to the supernova (SN) explosion that is associated with the 
formation of the NS, due to mass-tranfer phases. Any spin-orbit misalignement therefore is 
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expected to be introduced by the SN that forms the NS. iKalogeral (120001 ) argues that tilt 
angles greater than 30° are expected for 30% — 80% of coalescing BH-NS binaries, whereas 
tilt angles of 50° — 100° are expected for at most 70% of such systems. In order to account 
for these findings, we set up some simulations for misaligned mergers, covering a wide range 
of tilt angles from 30° - 180°. 



3.2.1. Setting Up Initial Conditions 

The selection of initial conditions for the equatorial mergers is straight forward: select 
an initial radius r -outside the Roche limit- and calculate the angular velocity needed 
for the relaxation procedure using Eqs. (1251) and (fTBl) . Obviously, Eq. (|25|) does not hold 
for non-equatorial orbits, since its derivation assumes that 9 = tt/2. When moving away 
from the equatorial plane, there is a third constant of the motion_| appearing (the Carter 
constant Q), as a result of the existence of a killing tensor. (The geodesies equations of the 
Kerr spacetime are included in Appendix B). In order to find a stable, so-called spherical, 
non-equatoria l orbit on which to initially place the NS, we follow the technique described in 
Hughesl ( 2000 ). (For an in depth analysis of the procedure, refer to Hughes 2000 and Wilkins 



1971). Here, we point out the basic steps fo r finding and s etting up numerically the initial 



conditions for such an orbit, as presented in iHughesI (120001 ) 



For a circular orbit R(= dr/dr) = = R', where the prime indicates differentiation 
with respect to r . Furthermore, for the orbit to be stable R" < should also hold. One 
can specify a unique orbit by fixing r and L z . The conditions R = = R' will fix the other 
two parameters of the orbit, E and Q. The inclination of the orbit (which is a constant of 
the motion) can be calculated as : 

cos(i) = ; Lz (26) 

(with i being zero on the equatorial plane where the Carter constant Q is also zero). 

The most bound orbit (in terms of energy) for a given radius r$ is the equatorial prograde 
orbit (in the same sense the retrograde orbit corresponds to the least-bound orbit). Therefore, 
one can start by choosing a radius r in the equatorial plane (where Q — 0) and then, by 
solving R — = R', calculate the energy E and angular momentum L z for that orbit. Solving 



the other two being the energy E and angular momentum L z , related to the stationary and axial Killing 
vectors respectively. 
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analytically the c ondition (with R being defined by Eq. (Bl)) one gets for the prograde and 



retrograde orbits (IHughesI |2000| ) 



E pro = )~ 2v2 +^ (27) 



y/l - 3v 3 + 2pv 3 
1 — 2pv 3 + p 2 v' 
y/l - 3v 2 + 2pv 3 
. — 2v 2 — pv 3 
1 — 3v 3 — 2pv 3 
1 + 2pv 3 + p 2 i 
\/l — 3v 2 — 2pv 3 



LP™ = - 1 " 2 ^ 3+PV (2t 



E ret = 1-2V 2 -PV 2 (29) 



L? = rv l+M±£t (30) 



where v = y/M/r and p = a/M. 



With E and L z numerically known for a given equatorial orbit of radius Tq, one can 
proceed into finding non-equatorial stable orbits. The way to do so is by keeping the radius ro 
fixed, decreasing the value of L z and solving again for the conditions R = = R', which now 
are going to give the energy and Carter constant. This way one can keep on decreasing the 
value of L z until it reaches the value L r z et or until R" = (marginally bound). The stability 
of the new, inclined orbit (with the inclination given by Eq. ff26l) ) should be checked with the 
requirement R" < 0. The angular velocity of this orbit can now be numerically determined 
as fl^ = d<fi/dt by using the geodesies equations for defi/dr and dt/dr (Eqs. (B3) and (B4)). 
The velocity for a stable, spherical, non-equatorial orbit of radius ro and inclination % will 
be v y = r Q^ cos(z) 

Since we have followed this numerical method to find initial conditions for the inclined 
mergers, tuning the r and L z values in order to solve for an exact value of inclination angle 
was not trivial to do. In any case we tried to get as close to the value of the desired inclination 
angle as it was numerically possible. For example, the inclined merger labelled to be 30° was 
in reality 29.6577°. In Table [2] we have listed the actual values of orbital inclinations used 
in our simulations, but when referring to them either in the text or in plot labels we use the 
rounded off values. 



3.2.2. Results 



We have run five simulations of non-equatorial mergers, for five different inclination 
angles all with a/M = 0.99 (Table [2]). The rest of the characteristics of the simulations 
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presented in this section, i.e., BH and NS masses, polytropic index T, number of SPH 
particles, are as mentioned in §3.1 . 

Runs 14 (i « 90°) and 15 (retrograde orbit) lead to the NS being entirely lost into the 
BH's horizon after spiraling around it for several orbits. The difference between the outcomes 
of those two mergers is that, for run 14 the feeding of the BH does not take place in the 
equatorial plane: the NS follows a 3-D spherical orbit of decreasing radius before it finally 
vanishes completely into the BH. 

The outcome of Run II is presented in Figs. [11] and [12j As in the equatorial mergers 
of high a/M presented in §3.1, the merger results in the formation of an expanding tail of 
the NS's material, with the innermost part of the helix 'feeding' the BH. For this case the 
helix does not lay on the equatorial plane , although the feeding point does. By the end 
of the simulation almost 40% of the NS mass remains outside the BH's horizon with most 
of it unbound and only about 6% bound. We followed the procedure described already in 
§3.1 and used the geodesic integrator to further evolve the results of the SPH simulation as 
a test of the validity of our definition of bound and unbound mass and to check the spatial 
distribution of the bound material around the BH. The whole unbound part of the surviving 
material escapes outwards, whereas the bound mass forms a stable torus outside the BH's 
horizon. Fig. [13] shows the spatial distribution of pericenters (in red - upper right and lower 
panels) and apocenters (in blue - upper left and middle panels) on the x-y and x-z planes. 
Every particle is depicted at the moment of its own apocenter (or pericenter) (i.e., those 
plots do not correspond to a snapshot of specific time of the equivalent merger). 

Run 12 corresponds to a merger with initial inclination of ~ 45°. As the infall of the 
NS's material into the BH starts, at about t/M = 1500 after the beginning of the simulation, 
an outwards expanding tail is forming (Figs. [T4] and [T51) and by the end of the simulation 
25% of the initial NS mass (Fig. [TB]) is left unbound and escaping. No bound, disk-forming 
material was resolved in this case. 

For an inclination of ~ 70° (Run 13) the outcome of the merger is no different than 
that of Runs 14 and 15: the whole NS disappears into the BH's horizon. The morphological 
features of this merger, however, are somewhere in between those of the low inclination 
mergers, where there is material surviving outside, and the higher inclination ones, which 
ended up being completely catastrophic for the NS. After the NS orbits around the BH 
following a spherical orbit of decreasing radius, it starts accreting into the BH with the 
feeding point being well above the equatorial plane. As the infall continues, the part of the 
NS that is still outside the BH's horizon moves towards the equatorial plane and at the last 
stages of the merger, when the feeding point has sufficiently approached the equatorial plane, 
a small, expanding, spiraling tail forms which is energetically unable to escape or expand 
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significantly and ends up being accreted into the BH as well. At the end of the simulation 
there is no surviving mass (Fig. [TTj) . 

The difference between equatorial and inclined mergers in terms of the surviving bound 
mass seems to be the vertical (z-axis) distribution of the NS debris. In the equatorial case 
the bound material remains in an equatorial disk (— 0.15M < z < 0.15M) whereas in the 
inclined case the disk is replaced by a torus with —1M < z < 1M, and this torus is also 
more massive than the disk in equatorial mergers. 

We note that the geodesic integrator runs are taken just as an approximation and a 
qualitative test of the further evolution of the two groups of surviving material (bound and 
unbound), for the cases where hydrodynamics is of no importance and the SPH particles can 
be treated as free particles; e.g., for unbound particles at great distances from the BH. 



4. Gravitational Wave Signals 

We present here the GWs extracted for some of the simulations described in the previous 
section. We have used the standard quadrupole formula, in which the two polarizations of 
the waveform are given by 

rh x = 2I xy (31) 

vh-\- — I X x lyy (32) 

where r is the distance to the observer, h x and h + the cross and plus polarization modes, 
and Iij the second time derivative of the second mass moment. 

Fig. [TS] shows the GWs extracted from runs El, E10, 12 and II. The upper left plot 
corresponds to Run El (equatorial merger of a maximally spinning BH). As the inspiral of the 
NS towards the BH continues, the amplitude and frequency of the waveform is increasing, to 
reach its maximum just before the NS gets disrupted and starts shedding mass into the BH. 
From this point on, the amplitude decreases very rapidly and the signal eventually vanishes. 
The GWs for an equatorial merger of a non-spinning BH (run E10 in Table 1) is shown in 
the upper right plot of Fig. [181 Again we see the chirp signal reaching its maximum just 
before the NS is lost into the BH's horizon. For a non-spinning BH the disruption of the NS 
is taking place very close to the horizon (as discussed above) and it just barely starts getting 
tidally distorted before plunging almost intact into the BH. As a result, we do not see the 
gradual diminishing of the wave amplitude as observed in the spinning case. The middle 
and lower panel plots of Fig. [TBI show the GWs for Runs 12 and II, two simulations for the 
merger of a maximally spinning BH with a NS set in an inclined orbit of inclination 45° and 
60°, respectively. Again we observe the characteristic chirp signal as the NS is inspiralling, 
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with the frequency of the signal decreasing as the wave amplitude is increasing and reaches 
its maximum just before the NS disruption takes place. Additionally, the precession of the 
orbital plane is seen as a secondary chirp signal of smaller amplitude. At this point, the 
main core of the newly disrupted NS is found orbiting on a different plane than the rest of 
the NS material that has just started forming an outwards moving tail. 

The GW signals shown in Fig. [18] are, at best, qualitatively correct, for two reasons. 
First, recall that the inspiral part of the calculations presented here is affected by the artificial 
strengthening of radiation reaction that we have used in order to drive the initial stages of 
our simulations (as described in §2.3). As a result, the early parts of the GW signals suffer 
a time compression that would not be present in reality. Second, the simple quadrupole 
formula of Eqs. ED and E21 applies strictly to the case of a quasi-Newtonian source, while in 
reality the NS fluid is moving highly relativistically in the strong field region around a BH. 
Our intention for now is merely to provide a qualitative preview of the GW signals produced 
by these BH-NS mergers, and not to derive quantitatively accurate waveforms. As part of 
our future work we plan to develop a more sophisticated and more accurate treatment of the 
radiation reaction and GW extraction, which will then allow us to study in greater detail 
how the merger's parameters produce distinctive features in the GW signals and energy 
spectrum. 



5. Discussion and Summary 

We have performed simulations of BH-NS mergers with mass ratio q ~ 0.1 and poly- 
tropic index T = 2 using a 3-D relativistic SPH code. We investigated equatorial mergers 
for various values of the BH spin and we also carried out simulations for inclined mergers 
(considering a tilted orbital plane for the NS with respect to the BH's equatorial plane) in 
the case of an extremal Kerr BH. We find that the outcome of the merger depends strongly 
on the spin of the BH and the inclination of the orbit. More specifically, for equatorial 
mergers, the survival of NS material is possible only for mergers with a/M > 0.7 and, in 
that case, the percentage of surviving material increases with increasing BH spin, varying 
from about 33% to ~ 1% for a/M decreasing from 0.99 to 0.7. Complete disruption of the 
NS happens for all values of a/M. Most of the surviving material gets ejected from the 
vicinity of the BH and only a few percent stays bound to the BH, forming a relatively thin 
stable disk outside the BH's horizon. Only for very high BH spins (a/M > 0.95) do we see 
a substantial fraction (~ 3%) of the disrupted NS material remaining bound. 

The outcome of inclined mergers (for fixed a/M = 0.99) shows strong dependence on the 
value of the inclination angle. For sufficiently low inclinations (< 60°) there is always a large 
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fraction of the NS mass surviving the merger, up to ~ 40%, depending on the inclination 
angle. Moreover, for these mergers the formation of a thick stable torus of substantial bound 
mass is also strongly inclination-dependent. Whenever the inclination exceeds 40 — 45° the 
fraction of bound material drops to levels unresolvable by our present calculations, although 
there is still a significant fraction (~ 25%) of the material that is being ejected. 

Given the sensitive dependence of our results on the BH spin, one might wonder whether 
the accretion of NS material onto the BH could lead to significant spin up, with a corre- 
sponding change in the spacetime metric which is not accounted for in our SPH code (where 
the BH spin parameter is assumed fixed during the entire evolution). For each one of the 
runs presented in §3 we have calculated the change in the BH spin by adding to the initial 
angular momentum of the BH the angular momentum of the total accreted mass (a poste- 
riori). For Runs E1-E5 the total change in the BH spin after the accretion was less than 
0.1. More precisely, for those runs with 0.7 < a/M < 0.95 the final BH was spun up by 
less than 0.1, while for Runs El and E2 the final BH spin decreased to a/M = 0.965 and 
a/M = 0.946, respectively. Those changes are small enough that we think they would not 
have affected the final results of the simulations. For the runs with BH spin < a/M < 0.6 
the change in the final BH spin varies between 0.27 for a non-spinning BH and ~ 0.1 for a 
BH with initial a/M = 0.6. Even for these cases the change in the BH spin will not alter 
the result of the merger, as all simulations with BH spin < 0.7 resulted in the accretion of 
the entire NS mass by the BH. 

Another possible caveat concerns our Newt onian treatment of t he NS fluid self-gravity. 



We should stress here the point made already in lFaber et al.l (l2006bl ). that a more accurate, 
relativistic treatment of the NS fluid's self-gravity could lead to even faster mass transfer 
to the BH, as it might result in the expansion of the NS on a dynamical scale during the 
mass loss phase. Although this effect could have significant impact on the outcome of some 
mergers, we have observed in all our simulations that the complete disruption of the NS lasts 
typically less than one dynamical time from the onset of mass loss and therefore we believe 
that the adoption of a simple Newtonian treatment of self-gravity in our code is justified. 

Our findings have direct implications for the viability of BH-NS mergers as progenitors 
of short GRBs. The possibility of short GRBs being powered by such mergers depends on 
the characteristics of the binary, mainly the BH spin and the orbital inclination. Although 
not all BH-NS mergers seem to be promising as progenitors for short GRBs, binaries with 
fast spinning BHs and zero to moderate inclinations appear to form a disk or torus around 
the BH, massive enough to power a short GRB when later accreted by the BH. Therefore, 
only a fraction of BH-NS binaries could be associated with the production of short GRBs. 
The question of whether this fraction is high enough to explain the observed short GRBs 
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has yet to be answered as it depends on uncertain quantities such as the BH spin and mass 
distributions, as well as the orbital inclination distribution for BH-NS binaries. Current 
population synthesis codes try to answer such questions and give better constraints on the 
distributions of t hose quantities. A recent discussion of this question in connection to short 



GRBs is given in iBelczynski et al.l ( 120071 ). The fact that the short GRB progenitors need 
not belong to a homogeneous group and that NS-NS binary mergers could be equally viable 
progenitors of short GRBs is also an open issue currently under investigation. 

The complete tidal disruption of the NS in all the merger calculations that we presented 
in this paper is in q ualitative agreement with previous Newtonia n studies of BH-NS merg- 



ers with soft EOS (IRosswog et al.l |2004| ; iLee &: Kluzniakl Il999bl ). In our simulations that 



involved a Schwarzschild (non-spinning) BH, no NS material was observed being ejected or 
forming a disk around the BH. The accretion of the whole NS onto the BH is very prompt 
for this case. This is to be expected, as for the particular mass ratio (q ~ 0.1) considered 
in our simulations the tidal limit is found inside the ISCO for a/M = and thus we would 
not expect to find a stable disk forming around the BH. The ejection of material and the 
formation of an accretion disk were observed only for mergers with highly spinning BHs. 

In Paper II we will extend our simulations to include different mass ratios, although still 
in the limit where the use of a fixed background is justified (i.e., for BH masses much larger 
than the NS mass), and we will explore the effects of changing the NS EOS and relaxing 
our assumption of an initially synchronized NS. In this context we would like to add a brief 
comment about the importance of the NS EOS for the dynamics of the merger right after 
the disruption of the NS. Fig. [T9] shows the averaged ratio of the pressure gradient forces 
over the gravity forces as a function of time (for the run with higher resolution presented 
in §3.1). We see clearly that right after the disruption of the NS (at t/M ~ 400) the 
pressure forces become insignificant and the decompressed NS material essentially follows 
ballistic trajectories, indicating that the details of the EOS no longer play any role. Thus 
one need not worry that the decompressed material should be described by a different EOS. 
Only the high-density EOS describing the interior of the NS prior to disruption may have a 
significant effect on the outcome of the merger. Ultimately, the adoption of a full numerical 
relativistic scheme, where no simplifying assumptions on the background spacetime metric 
are made, would be ideal, as it would allow us to explore mergers for binaries with higher mass 
ratios. Although the many Newtonian studies of these mergers have provided a qualitative 
overview of their possible outcomes, there is a consensus that relativistic effects will play 
a very significant role in the hydrodynamics of the NS disruption (as already shown by 



Shibata fc Uryul 12006k for the non-spinning BH case), and therefore in understanding the 



GW emission and the possible connection of such mergers to short GRBs. 
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A. K-S form of the Kerr metric 



Kerr presented his solution for the first time (IKerrlll963l ) in the following format 



ds 2 = (r 2 + a 2 cos 2 9) (d9 2 + sin 2 9d(p 2 ) + 2(du + a sin 2 9d(f>) x (dr + a sin 2 

2Mr 
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f 2 + a 2 cos 2 6 



x {du + a sin 



(Al) 



By using: 



u — t + r 
r cos 9 = z 
(f — ia)e 1 ^ sin 9 = x + iy 



(A2) 



Eq. (|A1| can be transformed to an asymptotically flat coordinate system (lKerrlll963f ) 
First let's apply u = t + r on Eq. (IA1I) . This will lead to the more familiar form 



Zj Zj 



where 



and r is defined by : 



(r 2 + a 2 + — — a 2 sm 2 6>) sin 2 6><i0 2 + -^-dt dr 

4Mf 2Mf 
H — — — a sin 2 9d(j) dt + (1 — — )2a sin 2 9d<p dr 



£ = f 2 + a 2 cos 2 

i;4 /„2 „2\~2 2 2 



~4 12 2\~2 2 2 n 

r — (p — a )r — a z = 
with p 2 = x 2 + y 2 = z 2 . 

Now applying the rest of the transformations in (A2), one gets 



(A3) 



(A4) 



where 



2Mr 3 

ds 2 = dx 2 + dy 2 + dz 2 - dt 2 + — ^{k) 2 

r 4 + a l z 2 



jixdx + ydy) + a(xdy — udx) z , 
k = \- ^4 5—^ + ~dz + dt]. 

yZ _|_ /■) Z <y> 



r l + a 



(A5) 



(A6) 
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Let us make a useful parenthesis here: to understand better why the transformation 
rules (A2) were chosen, it is constructive to start with Eq. flAlj) and bring it to the following 
form 

ds 2 = \-du 2 + 2du df + Y,d6 2 + 2a sin 2 Bdr d<p + (r a + a 2 ) sin 2 Odcf 2 ] 

2 Mr 

+ ^-L(du + asm 2 dd<j) 2 ). (A7) 

Eq. (1A7P can be interpreted as following: the terms not containing the mass m, give the 
flat space metric in some coordinate system, while the last term can be expressed in terms 
of the null (tangent with respect to r] a p) vector, l a 

—l a dx a = du + a sin 2 Odcf) 
and therefore the line element ds 2 can be written in the form 



2Mf 

ds 2 = (ds 2 )fi at + — (l a dx a ) 2 



and the metric is 



9af3 — Va/3 + 



2Mr 



l a lp. 



(A8) 
(A9) 



This is actually the original way which Kerr followed to discover his solution^ 



The idea now is to find those transformations that will take the part of Eq. ( 1A7j) that is 
contained in the brackets to the standard representation of the Minkowski space. This will 
lead us to the already mentioned transformation rules (A2). 

Notice that the metric (IA50 is analytical everywhere except at x 2 + y 2 = a 2 (or else at 
p = a and z = 0). 



B. Geodesies of the Kerr spacetime in B-L coordinates 



The Kerr geodesies equations (IMisner. Thorne fc Wheeler! 1 1973T ) 



S 2 ( ^\ = [E{r 2 + a 2 ) - aL z } 2 - A[r 2 + (L z - aE) 2 + Q] = R 



8 Any metric of the form g a p = rj a p + Hl a lp, with H a scalar and l a a null vector field, is called a 
Kerr-Schild metric. 



S 2 ( ^ V = Q - cot 2 6L 2 Z - a 2 cos 2 0(e - e 2 ) = 



E [ ^ ) = csc 2 ^ + aE i 7 -^- - ^ L ~- 



dr V A A 



2 _l „2<\2 

- a z sin^ 



(r + a ; 2 • 2 



A 



+ aL z [l- 



r 2 + a 2 
A 



C. for equatorial orbits in K-S coordinates 

The metric in K-S coordinates 

, 9 2Mf . , 9 „ „ 9 „ , ,~ AaMf sin 2 6> , ~, „ . 9 „ .„ ~. 
ds 2 = -(1 ^)c?w + Ec$ 2 + 2cfo dr dv d<p' - 2a sin 2 d<p 

,„ 2 2 2Mfa 2 sin 2 6> . 9 _ ~ /9 
+ (r 2 + a H ) sin 2 6 d(f) 2 



where 



~4 / 2 , 2 , 2 2\~2 22 n 

r — (x + y + z — a )r — a z = 



E = f 2 + a 2 cos 2 0. 
Restricting the metric on the 6 = n/2 plane 

, 9 . 2M X , 9 , , 4aM , ~, , ~. . 9 9 2Ma 2 x ~ /9 
(is 2 = -(1 —)dv 2 + 2dvdr —dvd<p' - 2adrd(f) + (r 2 + a H — )d(f) 2 . 

Two Killing vectors are associated with the v— and 0'— invariance or Eq. ( 1C4I) 

f= (i,o,o,o) 

7f= (0,0,0,1) 
and the defined conserved quantities 
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e = -i-u (C7) 

l = ff-u (C8) 

with e and / being the conserved energy per unit rest mass and conserved angular 
momentum per unit rest mass respectively, u is the four- velocity. 

From Eqs. (IC7I) and (IC8I) and using w r = 

e = -0 w ti« - <^t/ (C9) 

i = ^X + w«*- ( C1 °) 

From Eqs. (El) and flCTOl) 

„ cfc 1 r/ „ 2 2 2Ma 2 2Ma ;1 

M = 1Z = 5 I + a + ~ " ) e - -=-fl ( C11 ) 



c?r f 2 — 2Mr + a 2 f f 

_ 1 

dr f 2 — 2mr + a 2 



= -f = - . , J (l - — Z + — e]. (C12) 



is defined as 



^ " ^/rf r (f 2 + a 2 + 3M^) -m^l/e 1 J 



We need to substitute for 1/e in Eq. (1C13j) . Heading for that we can make use of the 
normalization condition 

u-u = -l, (C14) 

to find 

e 2 -l ldr M l 2 -a 2 (e 2 -l) Mil-aef, ,„ , 

Based on Eq. (10151) . one can define the effective potential 
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V eff = -— + 



M l 2 -a 2 (e 2 -l) M{l-aef 



r 2f 2 

For a circular orbit of radius f , the effective potential has a minimum at f 

(r ) = 



and also from Eq. (1C15j) 



Veff(r ) 



e 2 -l 



(C16) 



(C17) 



(C18) 



From Eqs. (1C17P and (1C18P one can solve for l/e: start from Eq. (1C17I) 



e// 



(e 2 - 1) + I 2 + ^(a 2 - a 2 e 2 + Z 2 ) 2 - 12(1 - ae) 2 M 2 

2M ' 



Now, from Eqs. (ICT8|) and dCT9|) solve for l/e 



(C19) 



Z ^ -f 3/2 (a 2 + f(f - 2M)) + aM 1 / 2 {a 2 + f(3f - 4M)) 
i ~ M x / 2 (a 2 M - f(f - 2M) 2 ' 

Substituting Eq. flC20l) in Eq. flCT3l) will give 

M 1//2 



;3 / 2 ± aM 1 / 2 



(C20) 



(C21) 



with the upper and lower signs corresponding to co-rotating and counter-rotating orbits. 
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Table 1. Runs and Results for Simulations of Equatorial Mergers 



Run 


El 


E2 


E3 


E4 


E5 


E6 


E7 


E8 


E9 


E10 


a/M 


0.99 


0.95 


0.9 


0.8 


0.75 


0.6 


0.5 


0.2 


0.1 





Total NS mass outside r + 


33 % 


32% 


26% 


4% 


1% 


0% 


0% 


0% 


0% 


0% 


Bound NS mass outside r + 


2.5% 


2% 


0% 


0% 


0% 


0% 


0% 


0% 


0% 


0% 



Table 2. Initial Conditions and Results for Inclined Mergers 



Run II 12 13 14 15 



a/M 0.99 0.99 0.99 0.99 0.99 

r (M) 12 11.5 12 13 10 

inclination angle (°) 29.6577 45.031 70.3989 89.9494 180 

Total NS mass outside r + 37% 25% 0% 0% 0% 

Bound NS mass outside r+ 6% 0% 0% 0% 0% 



1M 



R t =3R NS (M/M NS ) ,/3 : R NS =15km 

R t =2R NS (M/M NS ) 1/3 : R NS =10km 



10- 



10 
M/M 



100 



Fig. 1. — Tidal disruption limits for a 1.4M Q NS in circular orbit around a BH of mass M. 
The radius of the ISCO is also shown for the Schwarzschild (dashed line) and maximally 
co-rotating Kerr cases (dash-d otted line). The two solid lines brack et the tidal limit for NS 
with different spins and radii (jWiggins fc Laill2000l ; lLai et al.lll994l ). Orbits are assumed to 
be in the equatorial plane and prograde. 
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Fig. 2. — The future (solid blue line) and past (dashed red line) horizons of a BH for various 
values of the Kerr parameter a. The straight (dash-dotted green) line represents the ring 
(curvature) singularity of the K-S metric, which exists only in the equatorial plane. 
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Fig. 3. — A 'horizon-trapped' orbit for different values of a. The initial conditions for the 
particle are fixed and only a is varied in the four orbits. The red circles represent the two 
horizons. It is clear that the particle successfully crosses the outer horizon and is eventually 
trapped at the inner horizon, except for the a = M case where the two horizons coincide. 
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Fig. 4. — Test calculation for a WD orbiting a much more massive BH at r = 8 M. The mass 
ratio in this case is q ~ 4 x 10 5 . The constant radius of this circular orbit is maintained 
by our code to within better than 10~ 3 over one full period. Here each cross indicates the 
position of the WD (center of mass of all SPH particles) at a different time along the orbit 
(counter-clockwise) . 
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a/M = 



t/M = 9.97232 






(M) 



a/M = 



t/M = 60.021 1 




-5 
x (M) 
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t/M = 190.148 




-5 5 

x (M) 



a/K/ = 



t/M = 194.143 






x (M) 



a/K/ = 



t/M = 198.146 






x (M) 



Fig. 5. — Sequence of six snapshots from a simulation of a non-spinning BH-NS merger 
(equatorial projection). The 1.4 M NS gets completely disrupted by a non-rotating BH of 
mass M = 15 M & . The NS was initially placed outside the tidal limit (r ~ 8M for this NS 
of radius R ~ 15 km with a F = 2 polytropic EOS). The NS fluid is disappearing completely 
into the BH horizon (at r = 2M, indicated by the red circle) in a time t/M = 180 after the 
beginning of the simulation. 
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Fig. 6. — Same as in Fig. El but for a rotating BH with a/M = 0.99. Some of the NS material 
is still flowing inside the BH's horizon, but there is now also a tail of expanding material 
forming. By the end of the simulation (t/M ~ 550) the infall of material stops completely 
and about one third of the NS mass resides outside the horizon, with ~ 30% corresponding 
to unbound ejected material and the rest to the bound, disk-forming mass. 
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Fig. 7. — A close up view of the central part of four of the snapshots in Fig. [6] The three 
colored circles represent, in K-S coordinates, the BH's future horizon r + (blue), the past 
horizon r_(red) and the ring singularity at r = a (green). This illustrates what we discuss 
in §2.1.1: the metric's singularity at the inner (past) horizon leads to the illusory formation 
of a ring in the equatorial plane with radius r = r_. 
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Fig. 8. — Fraction of the NS's initial mass that resides outside the BH's future horizon as 
a function of time, for five different values of the BH's angular momentum a/M (0.75 < 
a/M < 0.99). The survival of NS material outside the BH horizon depends strongly on the 
BH's angular momentum. As a/M decreases the fraction of surviving mass decreases as well 
and the merger becomes completely catastrophic for the NS for values a/M < 0.7. Also note 
that the accretion of material through the BH horizon starts earlier as the BH's angular 
momentum decreases; this is to be expected as the horizon is moving further away from the 
BH with decreasing a/M. 
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Fig. 9. — Fraction of the NS's initial mass residing outside the BH's horizon as a function 
of time for five different simulations, corresponding to five different values of a/M (runs 
E1-E5). The black, blue and red lines correspond to total, bound and unbound NS material 
respectively. The top two plots correspond to a maximally spinning BH, using two different 
resolutions: 10 4 and 10 5 SPH particles for the left and right plots, respectively. For values 
of a/M < 0.95 the percentage of surviving bound material drops to unresolvable levels. 
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Fig. 10. — Spatial distribution of the apocenters (blue) and the pericenters (red) for the 
disk-forming SPH particles of the a = 0.99M equatorial merger (Run El). The upper left 
(right) and middle (bottom) panels correspond to the x-y and x-z plane projections for the 
apocenters (pericenters) respectively. The total mass of the disk is 2.5% of the initial NS 
mass. The mean pericenter value (inner radius of the disk) is r peri ~ 2.5M and the mean 
apocenter (outer radius of the disk) IS T apo ^ 

30M. 
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Fig. 11. — Snapshot of Run II (up) and its x-z projection (down) towards the end of the 
simulation. The red circle represents the BH's outer horizon. The material that survives the 
merger forms an expanding helix. For this particular run 30% of the initial NS mass gets 
dynamically ejected and 6% remains bound, forming a stable torus outside the horizon. 
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Fig. 12. — Fraction of the NS mass that resides outside the BH's future horizon for Run II. 
At time t/M = 1950 after the beginning of the simulation, the surviving material stabilizes 
at ~ 37%. 31% corresponds to unbound escaping mass, while 6% stays bound and forms a 
torus around the BH's horizon. 
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Fig. 13. — Spatial distribution of the apocenters (blue) and the pericenters (red) for the 
bound SPH particles of Run II. Note that in inclined mergers the bound mass forms a 
thick torus rather than a disk (as for equatorial mergers; compare Fig. [101) . The upper left 
(right) and middle (bottom) panels correspond to the x-y and x-z plane projections for the 
apocenters (pericenters) respectively. The total mass of the torus is 6% of the initial NS 
mass. The mean pericenter value (inner radius of the torus) is r per i ~ 2M and the mean 
apocenter (outer radius of the torus) is r apo ~ 25M. 
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Fig. 14. — Snapshot from Run 12 (upper panel) and its x-z projection (lower panel) towards 
the end of the simulation. The red circle represents the BH's outer horizon. By the end of 
the simulation an outwards expanding helix has formed, ejecting 25% of the initial NS mass 
away from the BH. The fraction of surviving bound material was too small to be resolved in 
this calculation. 
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Fig. 15. — Four 3-D snapshots from Run 12 towards the end of the simulation. The NS has 
already been disrupted and starts shedding material into the BH's horizon. An expanding 
helix of unbound material forms at the same time, resulting in the ejection of 25% of the NS 
mass. 
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Fig. 16. — Fraction of the NS mass that resides outside the BH's future horizon for Run 12. 
At time t/M = 1600 after the beginning of the simulation, the surviving material stabilizes 
at 25% almost all of which is unbound and therefore escaping outwards. 
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Fig. 17. — Fraction of the NS's mass that resides outside the BH's future horizon for Run 13. 
At tine t/M = 1640 after the beginning of the simulation, the whole NS mass has crossed 
the BH's outer horizon. 
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Fig. 18. — GW signals extracted from some of the simulations. The two upper plots cor- 
respond to runs El (left) and E10 (right). The middle panel plots are for run 12 with the 
right side showing a blowup of the left plot focussing on the later stages of the evolution. 
The lower panel shows the GW signal for run II. Again the right side shows a blowup of the 
plot on the left. 
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Fig. 19. — The averaged ratio of pressure gradient forces to gravitational forces as a function 
of time for our highest resolution run. Right after the disruption of the NS (t/M ~ 400) 
the pressure forces become negligible in comparison to the gravitational forces. As a result 
the fluid's EOS plays no role anymore in the dynamics of the merger and the decompressed 
material follows essentially ballistic trajectories. 



